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Abstract 

We provide here all the procedures in Mathematica which are needed for the computation 
of the analytic images of the strong coupling constant powers in Minkowski (%l u (s; nj) and 
2^ lob 0)) and Euclidean (A„(Q 2 ;n f ) and Af loh (Q 2 )) domains at arbitrary energy scales 
(s and Q 2 , correspondingly) for both schemes — with fixed number of active flavours 
rif = 3,4,5,6 and the global one with taking into account all heavy-quark thresholds. 
These singularity-free couplings are inevitable elements of Analytic Perturbation Theory 
(APT) in QCD [IH3J, and its generalization — Fractional APT jlHS], needed to apply the 
^ APT imperative for renormalization-group improved hadronic observables. 
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Program Summary 

Title of program: FAPT 
Available from: 

http : //theor . j inr . ru/~bakulev/f apt .mat/FAPT.m 

http : //theor . j inr . ru/~ bakulev/f apt .mat/FAPT_Interp .m 

Computer for which the program is designed and others on which it is operable: Any 
work-station or PC where Mathematica is running. 

Operating system or monitor under which the program has been tested: Windows 
XP, Mathematica (versions 5 and 7). 

No. of bytes in distributed program including test data etc.: 
47 kB (main module FAPT.m) and 4kB (interpolation module FAPT_Interp .m); 
21 kB (notebook FAPT_Interp.nb showing how to use the interpolation module); 
10 888 kB (interpolation data files: AcalGlobfi.dat and UcalGlobfi.dat with £ = 
1, 2, 3, 3P, and 4]Q 

Distribution format: ASCII 

Nature of physical problem: The values of analytic images A U (Q 2 ) and 2l„(s) of 
the QCD running coupling powers a"{Q 2 ) in Euclidean and Minkowski regions, 
correspondingly, are determined through the spectral representation in the QCD 
Analytic Perturbation Theory (APT). In the program FAPT we collect all relevant 
formulas and various procedures which allow for a convenient evaluation of A V (Q 2 ) 
and 2t„(s) using numerical integrations of the relevant spectral densities. 

Method of solution: FAPT uses Mathematica functions to calculate different spec- 
tral densities and then performs numerical integration of these spectral integrals to 
obtain analytic images of different objects. 

Restrictions on the complexity of the problem: It could be that for an unphysical 
choice of the input parameters the results are out of any meaning. 

Typical running time: For all operations the running time does not exceed a few 
seconds. Usually numerical integration is not fast, so that we advice to use arrays of 
precalculated data and apply then the routine Interpolate (as shown in supplied 
example of the program usage, namely in the notebook FAPT_Interp.nb). 



1 The notebook FAPT_Interp.nb and all interpolation data files are available from the same place in 
the form of the zipped archive FAPT_Interp.zip of the size 1844 kB. In order that Mathematica notebook 
FAPT_Interp.nb can use these precalculated data files one should place the directory .\sources\ with 
all data files in the same directory as the main file FAPT_Interp.nb. 
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1. Introduction 



QCD perturbation theory (PT) in the region of spacelike four-momentum transfer 
(Q 2 = —q 2 > — hereafter we call it the Euclidean region) is based on expansions 
in a series over the powers of effective charge (or running coupling constant) a s (Q 2 ), 
which in the one-loop approximation is given by ai 1} (Q 2 ) = (4ir/b )/L with b being the 
first coefficient of the QCD beta function, Eq. g-Q, L = ln(Q 2 /A 2 ), and A = A QC d 
is the QCD scale parameter. The one-loop solution okt\Q 2 ) has a pole singularity at 
L = called the Landau pole. The £-loop solution ai\Q 2 ) of the renormalization group 
equation ([2J has an £-root singularity of the type L~ 1 ^ at L = 0, which produces the 
pole as well in the £-order term d^a l s {Q 2 ). This prevents the application of perturbative 
QCD in the low- momentum spacelike regime, Q 2 ~ A 2 , with the effect that hadronic 
quantities, calculated at the partonic level in terms of a power-series expansion in the 
running coupling, are not everywhere well defined. 

Such a singularity appeared first in QED |8] and was named "ghost" due to the 
negative residue at the corresponding propagator pole. It was interpreted as an indica- 
tion that quantum field theory is self-contradictory. However, as was shown in [HI ITU] , 
it is only a hint about the PT inapplicability in the region where the expansion param- 
eter is not small. Appearance of such "ghost" singularities from a theoretical point of 
view contradicts the causality principle in quantum field theory [TOj HI] , since it makes 
the Kallen-Lehmann spectral representation impossible. It also complicates the deter- 
mination of the effective charge in the timelike region (q 2 > — hereafter we call it 
the Minkowski region). In a seminal paper by N. N. Bogoliubov et al. of 1959 [12] . the 
ghost-free effective coupling for QED has been constructed using the dispersion relation 
technique. 

After the very appearance of QCD many researchers tried to determine the QCD effec- 
tive charge in the Minkowski region, which is suitable for describing the processes of e + e~ 
annihilation into hadrons, as well as quarkonium and r-lepton decays into hadrons. Many 
such attempts used analytic continuation of the effective charge from the deep Euclidean 
region, in which perturbative QCD is known to work well, into a Minkowski one, where 
actual experiments were performed: a s {Q 2 ) —> a s (s = —Q 2 )- In 1982 Radyushkin [13] 
and Krasnikov and Pivovarov [H] using the dispersion technique of [12] suggested regular 
(for s > A 2 ) QCD running coupling in Minkowski region, the well-known 7r _1 arctan(7r/L). 

In 1995 Jones and Solovtsov using variational approach [15] constructed the effective 
couplings in Euclidean and Minkowski domains which appears to be finite for all Q 2 and s 
and satisfy analyticity integral conditions. Just in the same time Shirkov and Solovtsov [1] , 



using the dispersion approach of [12], discovered ghost-free coupling A\(Q 2 ), Eq. (25a) 



in Euclidean region and ghost-free coupling 2li(s), Eq. (25b), in Minkowski region, which 
satisfy analyticity integral conditions: 

At the one-loop approximation the last coupling coincides with the Radyushkin one for 
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s > A 2 . This way of making the QCD's effective charge analytic in the timelike region 
was rediscovered later within an approach of fermion bubble resummation by Beneke and 
Braun [TS], and also by Ball, Beneke, and Braun [T7J. Due to the absence of singularities 
in these couplings, Shirkov and Solovtsov suggested to use this systematic approach, called 
Analytic Perturbation Theory (APT), for all Q 2 and s. 

Recently the analytic and numerical methods, necessary to perform calculations in 
two- and three-loop approximations, were developed [T8H21"] . This approach was applied 
to the calculation of properties of a number of hadronic processes, including the width 
of inclusive r lepton decay to hadrons [25H2H], the scheme and renormalization-scale 
dependencies in the Bjorken [30], El] and Gross-Llewellyn Smith [32] sum rules, the width 
of T meson decay to hadrons [33J, etc. Moreover, APT was applied to the analysis of 
the processes with two scales rather than just a single scale, namely: the pion-photon 
transition form factor [34, 35] and the pion electromagnetic form factor in the 0(a s ) 
order [54H3B] . To summarize, we can say that APT (see reviews [3TH3U] ) yields a sensible 
description of hadronic quantities in QCD, though there are alternative approaches to the 
singularity of effective charge in QCD — in particular, with respect to the deep infrared 
region Q 2 < A 2 , where appearance of nonzero hadronic masses may be important |4"0Tl4"2"] . 
The main advantage of the APT analysis is much more faster convergence of the APT 
non-power series as compared with the standard PT power series, see in [13J HI]. 

Three-point functions, used in describing the pion electromagnetic form factor or 
7*7 — > it transition form factor, contain logarithmic contributions at the next-to-leading 
order of the QCD PT, related to the factorization scale. If one set the factorization scale 
proportional to the squared momentum-transfer, = Q 2 , then these logarithms will go 
to zero, but additional RG factors of the type (a s (Q 2 ) / a s (nl)) u , with v = 7 n /(2 ) being 
a fractional number, will appear in the Gegenbauer coefficients of the pion distribution 
amplitude. In both cases spectral densities, used to construct analytic images of hadronic 
amplitudes, should change. This observation led Karanikas and Stefanis [4"5l 14*6] to pro- 
pose the concept of analytization "as a whole" , meaning that one should construct analytic 
images not only of effective charge and its powers, but of the whole QCD amplitude under 
consideration. 

A QCD inspired generalization of APT to the fractional powers of effective charge, 
called Fractional Analytic Perturbation Theory (FAPT), was done in [U [6] (for a recent 
review see [17], for a recent generalization see [H]), followed by the application [5] to 
the analysis of the factorizable contribution to the pion electromagnetic form factor. The 
crucial advantage of FAPT in this case is that the perturbative results start to be less 
dependent on the factorization scale. This reminds the results, obtained with the APT, 
applied to the analysis of the pion form factor in the 0(a 2 ) approximation, where the 
results also almost cease to depend on the choice of the renormalization scheme and its 
scale (for a detailed review see [17] and references therein). The process of the Higgs 
boson decay into a bb pair of quarks was studied within a framework of FAPT in the 
Minkowski region at the one- loop level in [19] and at the three-loop level — in [6] . Results 
on the resummation of non-power- series expansions of the Adler function of a scalar, Ds, 
and a vector, Dy, correlators within FAPT were presented in [50J. The interplay between 
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higher orders of the perturbative QCD (pQCD) expansion and higher-twist contributions 
in the analysis of recent Jefferson Lab data on the lowest moment of the spin-dependent 
proton structure function, T^Q 2 ), was studied in jSI] using both standard QCD PT 
and (F)APT. FAPT technique was also applied to the analysis of the structure function 
F 2 (x) behavior at small values of x [52], [53]. All these successful applications of (F)APT 
necessitate to have a reliable mathematical tool for calculations of spectral densities and 
analytic couplings which are implemented in FAPTj^] 

In this paper we collect all relevant formulas which are necessary for the running of 
A V \L\ and 2t^[L] in the framework of APT and its fractional generalization, FAPT. We 
discuss their proper usage and provide easy-to-use Mathematica [56] procedures collected 
in the package FAPT. A few examples are given. Here we do not consider the inclusion 
of analytic images of logarithms multiplied by fractional powers of couplings, namely, 
[a s (Q 2 )] 1/ • [ln(<5 2 /A 2 )] m , which are needed for the full implementation of FAPT, — we 
postpone it to the next paper. 

The outline of the paper is as follows. In the next Section we present the main formulas 
of perturbative QCD which are needed for the running of the strong coupling constant up 
to the four-loop level. Section [3] contains the basic formulas of APT and FAPTj^] Finally, 
in Section [4] we describe the most important procedures of the package FAPT and provide 
an example of using this package to produce some numerical estimations. We hope that 
for most practical applications it should be sufficient. In the Appendix we supply the 
complete collection of the developed procedures. 

2. Basics of the QCD running coupling 

The running of the coupling constant of QCD, a s (/i 2 ) = a s [L] with L = ln[/t 2 /A 2 ], is 
defined through)^] 



2 This task has been partially realized for both APT and its massive generalization [?2] as the Maple 
package QCDMAPT in [SI] and as the Fortran package QCDMAPT_F in [S5J. Both these realizations are 
limited to the case of fixed number of active quarks Nf — 3 only, and use approximate expressions for 
the two- and higher-loop perturbative couplings, compare, for example, Eq. (33) in |54j and our Eq. (JTj). 

3 Note here that FAPT includes APT as a partial case for the integer values of indices. 

4 We use notations f(Q 2 ) and f[L] in order to specify the arguments we mean — squared momentum 
Q 2 or its logarithm L = ln(Q 2 /A 2 ), that is f[L] = /(A 2 • e L ) and A 2 is usually referred to rit = 3 region. 



da s [L] 
dL 



f3(a s [L];n f ) = - a s [L] ^2b k (n f ) 




(2) 
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where n/ is the number of active flavours. The coefficients are given by [5TH66 
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( is Riemann's zeta function, with values (2 = Tf 2 /Q and (3 ~ 1.202 057. It is convenient 
to introduce the following notations: 



b (n f ) 



, a(fi 2 ;n f ) = f3 f a s (fi 2 ;n f ) and c k (n f ) = 



6 fc (n/) 



Then Eq. (|2| in the Z-loop approximation can be rewritten in the following form: 



day) [L;nf] 
d~L 



l + ^2 c k(n f ) (a w [L;n/])' 



k>l 



(4) 



(5) 



In the one-loop (1 = 1) approximation (c&(n/) = 6^(71/) = for all k > 1) we have a 
solution 

a ( i)[L] = i (6) 

with the Landau pole singularity at L — > 0. In the two-loop (I = 2) approximation 
(cfc(n/) = 6fe(n/) = for all k > 2) the exact solution of Eq. ([2]) is also known [671 EH] 



a (2) [L;n/] 



" c i (™/J 



with 2^ [L] 



(7) 



1 + H/_!(^[L]) 

where W_i[z] is the appropriate branch of Lambert function. 

The three- and higher-loop solutions am[L;rif] can be expanded in powers of the 
two- loop one, 0(2) [L; n/], as has been suggested in [THl I22H2U 

a {£ )[L;n f ]=J2 C n ) ( a ml L , n f]T ■ 



n>l 



Coefficients cjp are known and can be evaluated recursively. We use in our routine for 
the three-loop coupling expansion up to the 9-th power included: 
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Figure 1: Left panel: Comparison of the standard three-loop coupling a^\Q 2 ) (solid blue line) with the 
three-loop Pade one a^ P \Q 2 ) (dashed red line). Right panel: Relative accuracy S(Q 2 ) — (ai 3P \Q 2 ) — 
a^\Q 2 ))/a^\Q 2 ) of the three-loop Pade coupling as compared with the standard three-loop one. 



As has been shown in [23] this expansion has a finite radius of convergence, which appears 
to be sufficiently large for all values of rif of practical interest. Note here that this method 
of expressing the higher-£-loop coupling in powers of the two-loop one is equivalent to the 
'tHooft scheme, where one put by hands all coefficients in /3-function, except bo and b\, 
equal to zero and effectively takes into account all higher coefficients 6j by redefining 
perturbative coefficients di (see for more detail in [69J). 

Another possibility for obtaining the "exact" three-loop solution is provided by the 
so-called Pade approximation scheme. It is based on the Pade-type modification of the 
three-loop beta function: 



/5(3P) Os) 

da(3P) [L] 



4n 



b + 



bi a s /(47r) 



dL 



—a 



(3P) 



[L] 



1 - b 2 as/iAnh) 
ci a( 3P )[L] 



1 + 



1 - c 2 a( 3P )[L]/ci 



(10a) 
(10b) 



The last equation can be solved exactly with the help of the same Lambert function (here 
the explicit dependence on nj is not shown for shortness): 



-cr 1 



C2 /4 + w^(4 p \l}) 



with 



(3P) 



[L] 



-l+C2/cf — L/ci 



(11) 



The relative accuracy o 
three-loop equation (5 
356 MeV) and better t 



this solution as compared with numerical solution of the standard 
with I = 3 is better than 1% for Q 2 > 2 GeV 2 (with A 3 3) = 
han 0.5% for Q 2 > 5 GeV 2 , cf. Fig.[TJ 



In the four-loop approximation we use the same Eq. <\8h with corresponding coefficients 



C (4) = C (3) + A (4) 



(12a) 
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Figure 2: Left panel: Comparison of the four-loop coupling ai (Q 2 ) (solid blue line) with the 
three-loop one ots (Q 2 ) (dashed violet line). Right panel: Relative accuracy <534(Q 2 ) = (ai 4 \Q 2 ) — 
a^\Q 2 ))/ai 4 \Q 2 ) of the three-loop coupling as compared with the four-loop one. 



and 



6 ' b 12 

A (4) = ~c\ c 3 _ 4 a c 2 c 3 11 c§ ( 4 ) _ cfc3 9 c 2 c 2 c 3 19 c 2 . c 3 _ 49 a c\ 
7 20 5 20 ' 8 ~ 30 20 r 3 120 ' 

A ( 4 ) = cfca _ 41 c\ c 2 c 3 _ 946 ci cjj c 3 134 c 2 c§ 149 c 2 cj 
9 ~ 42 140 315 35 504 1 J 

In the left panel of Fig. 2 we show both couplings, the four-loop ats(Q 2 ) (solid blue 

— ^ 

line), and the three-loop as (Q 2 ) (dashed violet line) with fixed number of active flavors 
rif = 4. We normalize both couplings to the same value a s (m 2 z ) = 0.119 at the Z-boson 
mass scale. Numerically, as can be seen in the right panel of Fig.|2j the relative deviation 
5 U (Q 2 ) = (af\Q 2 ) - af\Q 2 ))/af\Q 2 ) varies from 6% at Q 2 = 1 GeV 2 to 0.5% at 
Q 2 = 25 GeV 2 . We also compared the four-loop coupling, calculated in accord with 
Eq. Q, with coupling, calculated using package RunDec [70] with the same normalization 
a s (m 2 z ) = 0.119, — the relative deviation appears to vary from 0.2% at Q 2 = 1 GeV 2 to 
0.04% at Q 2 = 25 GeV 2 . 

2.1. Global scheme 

Here we consider the scheme of the so-called "global pQCD" in which the heavy-quark 
thresholds are taken into account. We follow here to Shirkov-Solovtsov approach [T| [T8| 120] 
with the following values of pole masses of c, b, and t quarks: m c = 1.65 GeV, = 
4.75 GeV and m t = 172.5 GeV. In the MS scheme of the standard pQCD one needs to 
match the running coupling values in Euclidean domain at Q 2 corresponding to these 
masses: M4 = m c , M5 = my, and Mq = m t . In order to implement these matching 
conditions we need to use the original QCD coupling 

^(Q 2 ;n f ) = — ^— a w (Q 2 ;n/) , (13) 
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where the indicator ^ signals about the loop order of the approximation we usej^] 
In what follows we use all logarithms L with respect to three- flavor scale A|: 

L(Q 2 ) = In (Q 2 /A 2 3 ) . (14) 

Recalculation to all other scales is realized with the help of finite additions: 

In (Q 2 /A 2 k ) = L(Q 2 ) + \ k with \ k = In (Af/A*) , (15) 

and Afc — the corresponding to the specified value Uf — k scale of QCD. We also define 
the corresponding logarithmic values at the thresholds M k (k — 4-r- 6): 



L fc (A 3 )=ln (M 2 /A 2 ) 



(16) 



All QCD scales Af, f — 4, 5, 6, we treat as functions of the single parameter, namely, the 
three-flavor scale A 3 : 



A f A/(A 3 ) with A 3 > A 4 (A 3 ) > A 5 (A 3 ) > A 6 (A 3 



(17) 



which should be defined from matching conditions for the running coupling at the heavy- 
quark thresholds. 

For an illustration we consider here the two-loop approximation with the running 
coupling af [L; n/] 



-4 7T 



with zw[L; n/] 



s [ ' /J & (n / )d(n / )[l + W_iOML;n/])] 
'l/ci(nj)) exp [—1 + ?7r — L/ci(n/)]. Then matching conditions are 

af [L 4 (A 3 );3] = af [L 4 (A 3 ) + A 4 ;4] ; (19a) 

(19b) 
(19c) 



af [L 5 (A 3 ) + A 4 ; 4] = af [L 5 (A 3 ) + A 5 ; 5] ; 



(2) 



a, 



(2) 



[L 6 (A 3 ) + A 5 ;5] 



a. 



(2) 



[L 6 (A 3 ) + A 6 ;6] . 



These relations define constants X k with = 4 6 as functions of variable A 3 , namely 



A fe ^Af (A 3 ), 

and, as a consequence, the continuous global effective QCD coupling 

«f b;(2) (Q 2 ,A 3 ) = af [L(Q 2 );3]e(Q 2 <M 2 4 ) 



(20) 



+ a 



(2) 



+ a 



(2) 



+ a 



(2) 



L(Q 2 ) + Af(A 3 );4] 0(M 4 2 <Q 2 <M 5 2 ) 
L(Q 2 ) + Af (A 3 );5] ^(M 2 <g 2 <M 6 2 ) 
L(Q 2 ) + Af (A 3 );6b(M 2 <g 2 ). 



(21) 



5 Note here that the dependence a^(Q 2 ; nj) on n/ is the consequence of Eq. (|5j), where for I > 1 one 
has n/-dependent coefficients Cfe(n/). 
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Here is the list of partial values of (A3), \\ (A3) and L/(A 3 ) with / = 4,5,6 for 
A 3 = 400 MeV: 



A? 


= 333 MeV , 


A(2) 
iv 5 


= 233 MeV , 


A(2) 
iv 6 


= 98 MeV ; 


(22a) 


Af 


= 0.367, 


x(2) 
A 5 


= 1.08, 


x(2) 
A 6 


= 2.82; 


(22b) 




= 2.197, 


U 


= 4.750, 


^6 


= 12.162. 


(22c) 



In our m-file we use the following realizations. The QCD scales are encoded as 
Al[A, rif], A2[A, rif], and A3 [A, n/] (in Mathematica capital Greek symbol A can be writ- 
ten as \ [CapitalLambda] ): 

\ [CapitalLambda] £[A, k] = A£[A, n f = k}= Af{A) , (£ = 1 4, 3P ; k = 4 6) , (23a) 

the threshold logarithms — as A£4[A], A£5[A], and A£6[A] (in Mathematica Greek symbol 
A can be written as \ [Lambda] ): 

\ [Lambda] £k[A] = X£k[A] = In (A 2 /A£[A, k] 2 ) , (£ = 1 4, 3P ; k = 4 6) , (23b) 

the running QCD couplings with fixed nf — as aBarlfQ 2 , rif, A], aBar2[Q 2 , rif, A], and 
aBa,i3[Q 2 , n f , A] (in Mathematica Greek symbol a can be written as \ [Alpha]): 

\[Alpha]Bar£[Q 2 ,n/,A] = «Bar£[Q 2 , n f , A] = a { s e) [\n(Q 2 / A 2 );n f ], (£ = 1 4, 3P) , (23c) 

and the global running QCD couplings as aGloblfQ 2 , A], aG\ob2[Q 2 , A], and 
«Glob3[g 2 ,A]: 

\ [Alpha] Glob£[Q 2 , A] = aGlob£[Q 2 , A] = af ob > {£) (Q 2 , A) , {£ = 1 4- 4, 3P) , (23d) 

To be more specific, we consider here an example. We assume that the two-loop a s 
is given at the Z-boson scale as a^fh^ml/A 2 ); 5] = 0.119. We want to evaluate the 
corresponding values of the QCD scales A 3 , A 4 , and A 5 and the coupling af^'^lQ 2 , A) 
at the scale Q = M5. We show a possible Mathematica realization of this task. 

In[l] := «FAPT.m; 

Comment: NumDef FAPT is a set of Mathematica rules in our package which assigns typical 
values to the physical parameters used in our procedures. 

In [2]:= {MZ = MZboson/ . NumDef FAPT , Mb=MQ5/ . NumDef FAPT} 
Out [2]= {91.19, 4.75} 

Comment: evaluation of L23= A^ from ai 2 \ln(m 2 z / A 2 ); 5] based on the explicit solution, 
Eq. g§, Eq.@. 

In [3] : = L23=lx/ . FindRoot [\ [Alpha] Glob2 [MZ~2 , lx] ==0 . 119 , {lx , . 1 , . 3}] 
Out [3]= 0.387282 
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(2) (2) (2) 

Comment: evaluation of L24= and L25= from L23= A3 based on Eq. (23a). 



In [4] : = {L24=\ [CapitalLambda] 2 [L23 , 4] , L25=\ [CapitalLambda] 2 [L23 , 5] } 
Out [4]= {0.321298, 0.224033} 

Comment: evaluation of af oh,<y2 \M^) from L23= Ag 2 ''. 

In [5]:= \ [Alpha] Glob2 [Mb"2,L23] 
Out [5]= 0.218894 

3. Basics of FAPT 

In the end of the previous section we used for the running QCD couplings with fixed 
rif the Bar notations — aBarl[Q 2 , rif, A], aBa,r2[Q 2 , n/, A], and aBax3[Q 2 ,rif, A]. We did 
it on purpose to have a direct connection to our previous papers on the subject |4T[6l |4T]. 
where we used the normalized coupling a(/i 2 ) = /?/ a s (/i 2 ), cf. Eq. Q. To be in line with 
these definitions, we also introduce analogous expressions for the fixed- Nf quantities with 
standard normalization, i.e., 

A(Q 2 ) = ^p, = (24) 

which correspond to the analytic couplings A v and in the Shirkov-Solovtsov terminol- 
ogy PP. 

The basic objects in the (F)APT approach are spectral densities pv\cr; nf) which enter 
the Kallen-Lehmann spectral representation for the analytic couplings: 

Jo <y + Q 2 J-oo 1 + exp(L - L a ) 

poo — {€} ( \ poo 

mL s ;n f }= / Pv [CX] nf) da = / pf\L c -n f \dL a , (25b) 

Js ° JLs 

It is convenient to use the following representation for spectral functions 

A[L > Uf] = n Im W [L - ^ ^ = ntf fR(e) [L;n f ])» ' (26) 
which is based on the module-phase representation of a complex number 

In the one-loop approximation the corresponding functions have the most simple form 



f { i)[L] = arccos f -^== j , R {1) [L] = 



7T 2 (2* 
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and do not depend on nj, whereas at the two- loop order they have a more complicated 
form 



R m [L;n f ] 
cp [2 )[L;n f ) 



ci(n/) 1 + Wi (z w [L — nr; n/]) 
-R { 2)[L;n f ] 



arccos 



Re 



(29a) 
(29b) 



1 + Wi (z w [L - Z7r;n/]) 

with W\ [z\ being the appropriate branch of Lambert function. In the three-loop approx 
imation we use either Eq. (|8]) and then obtain 

-l 



^>(3)[£] 



e i<P(2)[L] e ikip (2) [L] 



(3) 



arccos 



R { 3)[L] cos (<P(2)[L]) 
R(2)[L) 



+ J2 c {3) ^ ^ cos ^ V{2) ^ ' 



fc>3 



or Eq. (11) — and then obtain 

i-! + w.(4 P, [£-"]) 



^(3P) [L] 



Cl 



arccos 



Re 



l-Cca/^ + Wi^lL-iTr] 
In the four-loop approximation we use Eq. (]8| and then obtain 



(30a) 
(30b) 

(31a) 
(31b) 



i2(4) [L] 



, \- r (4) e <fc ^)W 



k>3 



arccos 



R(A)[L] COS (^(2)M) 



,(4) gffljjj cos(fcy? (2 )[L]) 

'fef* R\ 2) [L] 



(32a) 



(32b) 



Here we do not show explicitly the ri/ dependence of the corresponding quantities — it goes 



inside through R {2) [L] = R (2) [L; n f ], <p m [L] = m [L; n f ], Cf = C?>[n f ], C£> = Cf >/] 



t(3)r 



t(4) 



T (4)r 



Cfc = Cfc(n/) with = 1 -j-3, and z^ P ^[L] 



y (3P) 



[L; n/]. In the left panel of Fig. we show 



both spectral densities in comparison. On the right panel of this figure we show the 
relative deviation of the Pade spectral density from the standard one: one can see that it 
varies from +1% at L « —7, reduces to —2% at L ~ 0, and then reaches the maximum 
of +2% at L « 3.5. 

In accordance with Eq. (21) the global spectral densities are constructed through the 
ri/-fixed ones in the following manner: 

= [L a - 3] 9 (L a < L 4 (A 3 )) + pf \h a + X { 6 e) (A 3 ); el 9 (L 6 (A 3 ) < L a ) 



p^ oh [L a ,A 3 



+ P. 
+ P. 



La + Af (As); 4] #(L 4 (A 3 ) < L CT < L 5 (A 3 )) 
L a + A? (A 3 ); 5l 9 (L 5 (A 3 ) < L ff < L 6 (A 3 )) 



(33) 
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Figure 3: Left panel: Comparison of the standard three-loop spectral density p{ '[L] (solid blue line) 
with the three-loop Pade one p\ [L] (dashed red line). Right panel: Relative accuracy 5\ [L] = 
(p^ 3P ' ) [L] — [L] ) I p^ [L] of the three-loop Pade spectral density as compared with the standard three- 
loop one. 

with L a = \n(a/A^) and the corresponding global analytic couplings are 

Ai^[L,A 3 ] = / P^dL a , (34a) 

J-oc 1 + exp(L - L a ) 

/oo 
pi^° h [L a ,A 3 ]dL a , (34b) 

4. FAPT Procedures 

In our package FAPT.m we use the following realizations for the spectral densities. 
RhoBar£[L, rif, v\ returns £-loop spectral density p„ {£ = 1, 2, 3, 3P, 4) of fractional-power 
v at L = ln(C} 2 /A 2 ) and at fixed number of active quark flavors nf. 

RhoBar^L, k, u] = ffl [L;n f = k], (£ = 1 4 4, 3P ; k = 3 4- 6) , (35a) 

whereas RhoGlob£[L, u, A 3 ] returns the global £-loop spectral density slob [L;A 3 ] (£ = 

1,2, 3, 3P,4) of fractional-power u at L — \n(Q 2 /Al), cf. and with A 3 being the QCD 
rif = 3-scale: 

RhoGlob£[L, u, A3] = p^ );glob [L; A 3 ] , (£ = 1 4- 4, 3P) , (35b) 

Analogously, AcalBar£[L, nf, v\ returns £-loop [£ = 1,2, 3, 3P,4) analytic image of 
fractional-power v coupling Au [L;nj] in Euclidean domain, 

AcalBar£[L, k, v] = Af [L;n f = k], (£ = 1 4, 3P ; k = 3 6) , (36a) 

and UcalBar£[L, nf, v\ returns £-loop {£ = 1, 2, 3, 3P, 4) analytic image of fractional-power 
v coupling t&v[L,rif\ in Minkowski domain, 

UcalBar£[L, k, v] = [L; n f = k) , {£ = 1 4 4, 3P ; fc = 3 4 6) , (36b) 
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In global case AcalGlob£[L, A3] returns £-loop analytic image of fractional-power v 
coupling Au ^' glob [L, A 3 ] in Euclidean domain, 

AcalGlob^L, u, A3] = A^'' glob [L, A3] , (£ = 1 4- 4, 3P) , (37a) 

and UcalGlob£[L, i/, A3] returns £-loop analytic image of fractional-power v coupling 
21^ ' 9lob [L, A3] in Minkowski domain, 

UcalGlob£[L, 1/, A 3 ] = 2lj/ );9 ' ob [L, A3] , = 1 4- 4, 3P) . (37b) 

We consider here an example of using this quantities in case of Mathematica 7. We 
assume that the two-loop QCD scale A 3 is fixed at the value A 3 = 0.387 GeV defined at 
the end of section 12.11 

In[l] := «FAPT.m; 
In [2] := L23=0.387; 

We determine the value of the two-loop QCD scale L23APT = A^' APT in APT, corre- 
sponding to the same value 0.119 as before, but now for the global analytic coupling: 

In [3] : = MZ = MZboson / . NumDef FAPT 
Out [3]= 91.19 

In [4] : = L23APT=lx/ . FindRoot [AcalGlob2 [Log [MZ"2/lx"2] ,l,lx] 

== 0. 119, {lx, 0.35, 0.45}] ; 

Out [4]= 0.379788 

Now we evaluate the value of ^ 2);glob [L, L23APT] for L = -5.0, -3.0, -1.0, 1.0, 3.0, 5.0 
with indication of the needed time: 



In [5] : = 


{L0=-5. 


, AcalGlob2 [L0 , 1 , L23APT] }//Timing 


Out [5] = 


{0.734, 


{-5. , 0.929485}} 


In [6] : = 


{L0=-3. 


, AcalGlob2 [L0 , 1 , L23APT] }//Timing 


Out [6] = 


{0.421, 


{-3. ,0.786904}} 


In [7] := 


{L0=-1. 


, AcalGlob2 [L0 , 1 , L23APT] }//Timing 


Out [7] = 


{0.422, 


{-1. ,0.60986}} 


In [8] : = 


{L0=1. , 


AcalGlob2 [L0 , 1 , L23APT] }//Timing 


Out [8] = 


{0.437, 


{1. ,0.434041}} 


In [9] : = 


{L0=3. , 


AcalGlob2 [L0 , 1 , L23APT] }//Timing 


Out [9] = 


{0.469, 


{3. ,0.301442}} 
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-2 

Figure 4: Left panel: Graphics produced in Out[10] for _4[? ;glob [L, L23APT] as a function of L. Right 
panel: Graphics produced in Out[ll] for a^ 2);glob [L, L23APT] as a function of L. 



In [10] : = {L0=5 . , AcalGlob2 [L0 , 1 , L23APT] }//Timing 
Out [10] = {0 . 531 , {5 . , . 219137}} 

Now we create a two-dimensional plot of „4 2);glob [L, L23APT] and 2ll 2);glob [L, L23APT] for 
L G [—3, 11] with indication of the needed time: 

In[ll] := Plot [AcalGlob2[L, 1 ,L23APT] ,{L,-3, 11}, MaxRecursion->l] //Timing 
Out [11]= -[19.843, Graphics (see in the left panel of Fig.\,4)} 

In [12] := Plot[UcalGlob2[L,l,L23APT] ,{L,-3, 11} ,MaxRecursion->l] //Timing 
Out [12]= {14.656, Graphics (see in the right panel of Fig.\,4)} 



5. Interpolation 

The calculation of the spectral integrals is a computational task requiring a long time, 
especially if one is using the result in another numerical integration procedure. Therefore, 
it seems reasonable to pre-compute analytic images of couplings for a fixed set of argument 
values, consisting of N points for each argument. For example we will consider in what 
follows the case of A^' 9lob [L, u, Ag ]. We will be interested in the following ranges of 
arguments: L G [-5,5], A^ G [0.2,0.5], and v =G [0.5,1.5]. Then 

Lmin = —5 ; Lmax = 5 ; DL = (Lmax — Lmin) / (N — 1) ; 

z/min = 0.5 ; z/max = 1.5; Dz/ = (z/max — umin)/(N — 1) ; (38) 

Amin = 0.2 ; Amax = 0.5 ; DA = (Amax - Amin) /(N - 1) . 

The table of calculated values is generated by Mathematica using the following command 

DATA = Flatten[Table[{{Li, v h A k }, AcalGlobl[L 4 , v h A*]}, {i, N}, {j, N}, {k, N}), 2] 

where Lj = Lmin + (i — 1)DL, Vj = z/min + (j — l)Dz/, and A^ = Amin + (k — 1)DA. Then 
we save all calculated results in the file "AcalGlobli.dat": 
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In [2] : XY = N[DATA] ; {outFile = OpenWrite ["AcalGlobli . dat"] , 

Write [outFile, XY] , Close [outFile] } 
Out [2]: {OutputStream ["AcalGlobli. dat 11 , 15], Null, "AcalGlobli .dat"} 

After that we can read them and use interpolation to reproduce function A^P' 9lob [L, is, A^] 
in the considered ranges of arguments values: 

In [3] : DATA = Read ["AcalGlobli . dat"] ; 

AcalGlobllnterp = Interpolation [DATA] 

in order to select the appropriate value of N. Now we can analyze the accuracy of inter- 
polation. In Fig [5] we show the dependencies of interpolation errors on the number of 
the used points N. One can see, that using the interpolation at N = 6 for Av' 9lob and 
<y\?hal°b p r0 vided accuracy not worse 0.005%. 

In the previous case we investigated the dependence of the accuracy of interpolation 
on the number of points at fixed L, v and A 3 . Let us now consider how the accuracy of 
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Figure 5: Relative errors of the interpolation procedure for Au' glob (left panel) and 2lJ^^ 9iob (right 
panel) calculated at various loop orders with fixed L = 3.5, v = 1.1 and A3 = 0.36 GeV. 
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Figure 6: Relative error of the interpolation procedure for A 9 V 1 X 1 (left panel) and St^^J 1 (right panel), 
calculated at various loop orders with A 3 = 0.36 GeV for N = 11 number of points. 
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the interpolation depends on the L. These results are shown in Fig. |6j From the last 
figure one can see that the maximum error of interpolation corresponds to the region 
L = -r- 5. The error in A^Iq 1 ^ is less than in 2l£zf l ° b . In any case, using N = 11 points 
for interpolation of pre-computed data for each parameter L, v and A 3 provides an error 
less than 0.01 %. 

To obtain the results much faster one can use module FAPT_Interp .m which consists of 
procedures AcalGlob£i[L, u, A 3 ] and UcalGlob£i[L, u, A 3 ]. They are based on interpolation 
using the basis of the precalculated data in the ranges L = [—5; 13]; z/ 1_1o °p = [0.5; 4.0] and 
A^°p = [0.150; 0.300]; zA loop = [0.5; 5.0] and A^L° 3 P = [0.300; 0.450]; z/ 3 " 1oo p = [0.5; 6.0] 

and A^L°3 P = [0.300; 0.450]; zA loop = [0.5; 7.0] and A^° 3 P = [0.300; 0.450]. For example, 
in the four-loop case module FAPT_Interp .m contains procedures 

AcalGlob4i = Interpolation [Read [". \\sources\\AcalGlob4i . dat M ] ] ; 
UcalGlob4i = Interpolation [Read [". \\sources\\UcalGlob4i . dat M ] ] ; 

which should be used with the same arguments L, v, and A 3 as the original procedures 
AcalGlob£[L, u, A 3 ] and UcalGlob£[L, u, A 3 ]. They provide much faster results of calcula- 
tions with high enough accuracy: 



In[l] : = 
Out[l] = 


Timing [AcalGlob4i [1 , 
{0. , 0.39298} 


1.1, 


0.36]] 


In [2] : = 
Out [2] = 


Timing [AcalGlob4 [1 , 
{0.405, 0.392964} 


1.1, 


0.36]] 


In [3] : = 
Out [3] = 


Timing [UcalGlob4i [1 , 
{0. , 0.375421} 


1.1, 


0.36]] 


In [4] := 
Out [4] = 


Timing [UcalGlob4 [1 , 
{0.359, 0.375372} 


1.1, 


0.36]] 
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Appendix A. Numerical parameters 

Here we shortly describe numerical parameters used in the package. 



16 



First, in FAPT .m we use the pole masses of heavy quarks and Z-boson, collected in the 
set NumDef FAPT: 



MQ4 : 
MQ6 : 



1.65 GeV, 
172.5 GeV. 



MQ5 : 
MZboson 



M b -- 
M z 



4.75 GeV; 
= 91.19 GeV. 



(A.l) 



Note here that all mass variables and parameters are measured in GeVs. That means, 
for example, that in all procedures of our package the following value MQ4 = 1.65 is used. 
The package RunDec of [70J is using the set NumDef with slightly different values of these 
parameters (M c = 1.6 GeV, M b = 4.7 GeV, M t = 175 GeV, M z = 91.18 GeV). 

Second, we collect in the set setbetaFAPT the following rules of substitutions hi — > 
k{n f ), cf. Eq.g, 

2 

--n }) bl : b x 



bO 


bo 


— >■ 


b2 


b 2 


-> 


b3 


bs 


-)■ 



38 

102 



2857 5033 



2 

149753 



325 2 

tit H n f 

1 54 f 



(A.2) 



1078361 



+ 



6 

3564 



162 
6508 



n f + 



50065 



162 



27 



71/ + 



6472 



n. 



C[3]. 



1093 



tit 

729 f 



Here we follow the same substitution strategy as in [70], but our 6j differ from theirs bf KS 
by factors 4* +1 : b; L = 4* +1 bf KS . In parallel, the set setbetaFAPT4Pi defines substitutions 
which are more appropriate to determine coefficients Q(n/). 



h bi(n f )/(An 



Appendix B. Description of the main procedures 

Here we shortly describe the main procedures of our package which can be useful for 
practical calculations. 

• RhoBar£[L,Nf ,Nu]: 

general: it computes the £-loop spectral density p^[L a ,nf,v]; 

input: the logarithmic argument L=L a = ln[cr/A 2 ], the number of active flavors 
Nf=n/, and the power index Nu=z/; 

output: pW; 

example: In order to compute the value of the four-loop spectral density 
p (4) [3.95,4, 1.62] = 0.0247209 one has to use the command 
RhoBar4 [3 . 95 , 4, 1.62]. 

• RhoGlob£ [L , Nu , Lam] : 

general: it computes the £-\oop global spectral density p^' sloh [L a , u, A n/=3 ]; 
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input: the logarithmic argument L=L a = ln[a/A^ 3 ], the power index Nu=z/, and 
the QCD scale parameter Lam=A n/=3 (in GeV); 

output: pW;s lob ; 

example: In order to compute the value of the four-loop spectral density 

p (4); g iob[ 3>95) 1.62,0.350] = 0.0221662 one has to use the command 
RhoGlob4[3.95, 1.62, 0.35]. 
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• AcalBar£[L,Nf ,Nu]: 

general: it computes the £-loop n/-fixed analytic coupling Av\L,nf] in Euclidean 
domain; 

input: the logarithmic argument L=\n[Q 2 / A 2 ], the number of active flavors Nf=n/, 
and the power index Nu=z/; 

output: A^P; 

example: In order to compute the value of the three-loop spectral density 
.4S 2 [3.95, 4] = 0.11352 one has to use the command 
AcalBar3[3.95, 4, 1.62]. 

• UcalBar£[L,Nf ,Nu]: 

general: it computes the £-loop n/-fixed analytic coupling 2ii^[L, rif] in Minkowski 
domain; 

input: the logarithmic argument L=ln[s/A 2 ], the number of active flavors Nf=nj, 
and the power index Nu=z/; 

output: 2t^; 

example: In order to compute the value of the three-loop spectral density 
9lS 2 [3.95, 4] = 0.1011 one has to use the command 
UcalBar3[3.95, 4, 1.62]. 

• AcalGlob£[L,Nu,Lam]: 

general: it computes the £-loop global analytic coupling A^' slob [L, v, A n/=3 ] in 
Euclidean domain; 

input: the logarithmic argument L=L a = ln[<r/A 2 /=3 ], the power index Nu=z/, and 
the QCD scale parameter Lam=A n/=3 (in GeV); 

output: Al £y ' gloh ; 

example: In order to compute the value of the two-loop analytic coupling 
^L6? 0b [3-95, 0.350] = 0.103858 one has to use the command 
AcalGlob2[3.95, 1.62, 0.35]. 

• UcalGlob£[L,Nu,Lam]: 

general: it computes the £-\oop global analytic coupling 2l^' gIob [L, u, A n/=3 ] in 
Minkowski domain; 

input: the logarithmic argument L=ln[s/A 2 /=3 ], the power index Nu=z/, and the 
QCD scale parameter Lam=A n/=3 (in GeV); 

output: 2l!f );glob ; 
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example: In order to compute the value of the two-loop analytic coupling 
2lS 2 ^ glob [3.95, 0.350] = 0.0932096 one has to use the command 
UcalGlob2[3.95, 1.62, 0.35]. 

All A n/=3 are in GeV, all squared momentum transfer Q 2 (Euclidean), central-of-mass 
energy squared s (Minkowski), and spectral-integration variables a are in GeV 2 . The 
number of loops I is everywhere specified in the name of a procedure. 
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